Electrical response of molecular systems: the power of self-interaction corrected 
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The accurate prediction of electronic response properties of extended molecular systems has been 
a challenge for conventional, explicit density functionals. We demonstrate that a self-interaction 
correction implemented rigorously within Kohn-Sham theory via the Optimized Effective Potential 
(OEP) yields polarizabilities close to the ones from highly accurate wavefunction-based calculations 
and exceeding the quality of exact-exchange-OEP. The orbital structure obtained with the OEP-SIC 
functional and approximations to it are discussed. 
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Gaining microscopic insight into the quantum- 
mechanical electronic effects that govern energy- and 
charge transfer in processes like light-harvesting, charge- 
separation in organic solar cells, or the response of 
molecular opto-electronic devices would be extremely 
beneficial to the understanding of these phenomena. 
But the computational complexity of solving the many- 
electron Schrodinger equation leaves little hope that 
wave-function based approaches can address these prob- 
lems any time soon. The formulation of quantum me- 
chanics without wavefunction, i.e., density functional 
theory (DFT) in the Kohn-Sham framework, is compu- 
tationally much more efficient and allows to handle sys- 
tems with up to several hundreds of electrons. Therefore, 
it appears as the ideal tool for investigating the above 
mentioned problems. However, the predictive power of 
DFT calculations depends crucially on the approxima- 
tions made in the description of the exchange-correlation 
effects. Structural, ground-state molecular properties 
are obtained with reasonable to excellent accuracy us- 
ing standard, explicit density functionals like the Lo- 
cal Spin Density Approximation (LSDA) or Generalized 
Gradient Approximations (GGAs). But these function- 
als notoriously fail in the description of charge transfer 
processes [l|, B and associated_problems like predicting 
the response [3j or transport (J| properties of extended 
molecular systems. There is, thus, a serious need for 
exchange-correlation approximations that allow to calcu- 
late response properties like polarizabilities of extended 
systems reliably on a quantitative scale and with bear- 
able computational costs. 

It has been demonstrated that improvements in the 
density-functional description of the response of conju- 
gated polymers can be achieved based on current den- 
sity functional theory H and related ideas Q, or by in- 
corporating full [!, 0, 11] or partial exact exchange. 
It has also been argued that correlation effects play a 
non-negligible role in the proper description of response 
properties [13]. However, evaluating the Fock integrals 
in exact exchange approaches increases numerical costs 
substantially, and the computational complexity of ap- 



proaches using exact exchange with "compatible" corre- 



lation is significant 11 1 



In this manuscript we demonstrate that these prob- 
lems can be overcome with a self-interaction correction 
(SIC) employed rigorously within Kohn-Sham theory. In 
the SIC-scheme, only direct, i.e., self-exchange integrals 
need to be evaluated, thus computational costs are low- 
ered. OEP-SIC yields highly accurate results for the re- 
sponse of extended molecular systems without involving 
empirical parameters. 

The first "modern" SIC was proposed by Perdew and 
Zunger as a correction to LSDA 13]. They devised the 
LSDA-SIC functional 
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where ££ C SDA is the LSDA exchange-correlation energy 
functional, £h the Hartree (classical Coulomb) energy, 



and m 



I the up- and down-spin densities, respectively, 
Aj and N± the numbers of occupied spin-orbitals, and 
rij )CT the orbital spin densities. Eq. ([1]) is not the only 
way in which a SIC can be defined [14j . but it is plau- 
sible and straightforward: The spurious self-interaction 
effects that are contained in the Hartree energy and the 
LSDA functional are subtracted on an orbital-by-orbital 
basis. However, a subtlety is buried in this seemingly 
simple equation: The functional depends on the orbitals 
explicitly, i.e., it is not an explicit density functional. The 
traditional way of approaching this problem has been to 
minimize the total energy with respect to the orbitals 
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Hohenberg-Kohn theorem, but it is outside the founda- 
tions of Kohn-Sham theory: minimizing with respect to 
the orbitals leads to single-particle equations with or- 
bital specific potentials instead of a global, local Kohn- 
Sham potential for all orbitals. But the existence of 
a common, local potential is one of the features that 
makes Kohn-Sham DFT attractive: Only with a local 
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potential is the non-interacting kinetic energy density a 
well-denned density functional, a local potential consid- 
erably simplifies numerical efforts, it facilitates the inter- 
pretation of results, and it yields not only corrected oc- 
cupied eigenvalues, but also corrected unoccupied ones. 
But on the other hand, Perdew-Zunger SIC [13|] has be- 
come popular in some areas of solid state physics ex- 
actly for the reason that it does not work with one com- 
mon local potential but with several orbital specific ones, 
because orbital specific potentials straightforwardly al- 
low to take into account orbital localization effects: SIC 
with orbital-specific potentials can treat, e.g., p- and d- 
orbitals of a crystalline solid on a different footing. In 
this way, Perdew-Zunger SIC can naturally distinguish 
between localized and delocalized states. In order to ben- 
efit from the advantages of working with a local poten- 
tial without loosing the ability to describe localization 
effects, schemes have been devised which make use of the 
fact that Eq. JI) is not invariant under transformations 
of the orbitals that change the individual orbital densi- 
ties but leave the total density unchanged. Calculating 
orbitals from a common Hamiltonian and then subject 
these orbitals to localizing transformations has proved to 
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be a successful scheme for solids 

However, localizing orbital transformations can be- 
come computationally involved in large finite systems 
and in time-dependent calculations. Therefore, yet an- 
other variant of the SIC has become popular. It uses 
the Krieger-Li-Iafrate (KLI) construction [12] to obtain 
the KLI-potential corresponding to Eq. {l} and evaluates 
Eq. J!) directly with the KLI orbitals [H El, [3, 
This approach has been justified as an approximation to 
the OEP-version of SIC (OEP-SIC), which is defined by 
evaluating Eq. ([T) with the orbitals obtained from the 
SIC Kohn-Sham potential that follows from the Opti- 
mized Effective Potential (OEP) formalism ll|. But to 
the best of our knowledge, the performance of the OEP- 
SIC approach itself has remained largely unexplored, and 
tests of the KLI- SIC approach were restricted to spheri- 
cal atoms [l7j . In this manuscript we demonstrate that 
OEP-SIC, but not KLI-SIC, allows to predict electric re- 
sponse coefficients of molecular systems very reliably and 
may thus become an important tool to investigate charge- 
transfer questions. 

One of the most demanding tests of a functional's abil- 
ity to adequately describe charge transfer is calculating 
the polarizability of hydrogen chains. It has been shown 
that obtaining the response of hydrogen chains correctly 
is even more difficult than obtaining the response of real 
polymers like polyacetylene Therefore, calculating 
the polarizability of hydrogen chains has become a bench- 
mark test for many-particle approaches from both the 
density functional [3, B Sj 0, an d the wave-function 
worlds [H, [29} . Since a response quantity like the po- 
larizability determines how a system reacts to a field 



that induces a density shift, calculating the polarizabil- 
ity also probes the ability to correctly describe charge 
transfer. As a second, positive side effect, investigating 
hydrogen chains also allows us to address the question of 
size-consistency of the OEP-SIC functional [11 El, Ell- 

Our calculations are based on a real space approach 
[30I ] which we employed to calculate the ground-state of 
hydrogen chains with alternating interatomic distances 
of 2 and 3 ao using KLI-SIC. From the converged KLI- 
SIC solution we calculated the true OEP following the 
iterative procedure described in [3ll ] , which is more cum- 
bersone for the SIC-LDA functional than for pure ex- 
change, but does converge. The ground-state calcula- 
tions (no electrical field applied yet) lead to a remark- 
able result. For the sake of clarity we discuss it using 
the specific example of the shortest chain, H4. The KLI- 
solution is spatially symmetric as expected and as de- 
picted in the left part of Fig. [TJ It is also manifestly 
spin-unpolarized, i.e., the self-consistent KLI iteration re- 
turns to a spin-unpolarized solution from a spin-polarized 
starting guess. But starting from the spin-unpolarized 
KLI solution and iterating the OEP to self-consistency 
without restriction on the spin polarization, we observe 
a spontaneous change in symmetry. After a few iterations 
of the OEP self-consistency cycle, the up- and down-spin 
orbitals separate and each orbital starts to center around 
one nucleus. For the interatomic distances of 2 and 3 ao 
frequently used in the literature 0, H, [28, 32 1, the effect 
is moderate but clearly visible, as shown in the middle 
of Fig. |T). If the interatomic distances are increased 
further, e.g., to 2.5 and 3 ao as shown in the right part 
of Fig. |T), the orbital localization becomes pronounced 
and one can undoubtedly associate one orbital with one 
nucleus. This effect is not only observed for H4, but also 
for the other hydrogen chains we studied. 

A conclusion from this finding is that the KLI-SIC po- 
tential is not necessarily a good approximation to the 
OEP-SIC potential. In order to understand this one 
should recall that the KLI-potential is justified as a mean- 
field approximation [ll, H HJ: The difference between 
the true OEP and the KLI-potential is a term of the kind 
(l/n(r))V • f(r), where f(r) is a well defined function 
depending on the full spectrum of Kohn-Sham orbitals. 
Averaged over the density, the term vanishes [2, El- 
But implicitly this mean-field argument assumes that the 
"averaged" term has little influence in the self-consistent 
iteration so that the density obtained with and without 
the neglected term are very similar. However, our calcu- 
lations show that this is not the case for the SIC func- 
tional: taking into account the term that is neglected in 
the KLI potential drives the self-consistent iteration to 
a very different solution. This is possible because the 
neglected term contains all orbitals and is thus relevant 
for unitary (in)variance and greater variational freedom. 
The breakdown of the KLI-SIC approximation may be 
a surprise in view of its good performace for atoms 17j . 




FIG. 1: Left: Orbital densities of H4 with interatomic distances of 2 and 3 bohr (ao), respectively, obtained from self-consistent 
KLI-SIC calculation. Up- and down-spin orbitals are identical. Middle: Spin-orbital densities for the same system obtained 
from self-consistent OEP-SIC calculation. Right: Spin-orbital densities of H4 with interatomic distances of 2.5 and 3 ao, 
respectively, obtained from self-consistent OEP-SIC calculation. The orbital localization increases with increasing interatomic 
distance. 



but appears less surprising in view of other drawbacks 

The hydrogen chain ground-state results also naturally 
trigger thoughts about the bulk limit that one would 
reach by adding ever more atoms. We briefly want to 
ponder this case. Recall that for an infinite lattice of hy- 
drogen atoms with a lattice constant that tends to infin- 
ity, the exact Kohn-Sham orbitals are delocalized Bloch 
orbitals for which the self-interaction correction vanishes 
on a per-atom basis [l3j]. Using such orbitals in Eq. (JTJ) 
yields the (wrong) uncorrected LSDA energy. Inherent 
to the logic of this argument is a certain order of taking 
the two "infinity limits" : first the number of atoms tends 
to infinity, i.e., an infinite lattice is built, and then the 
lattice constant is made ever larger. 

Our calculations suggest that a different result is ob- 
tained if the order of taking these two limits is inter- 
changed. For finite systems of largely separated hydro- 
gen atoms, our OEP-SIC calculations lead to localized or- 
bitals and thus, a self-interaction corrected energy. Now 
imagine building up an ever larger lattice of hydrogen 
atoms with ever larger interatomic separation by adding 
atoms to a finite starting system. At each step of this 
buildup process, one will be dealing with a large but fi- 
nite system. Our calculations suggest that at each stage 
of the build-up process, OEP-SIC will yield localized or- 
bitals and thus a self-interaction corrected energy. This 
idea is in line with earlier findings that revealed that it 
makes a great difference whether the surface of an ex- 
tended system is explicitly taken into acoount or not 
[3^ |. In any case our results show that OEP-SIC can 
yield localized orbitals that differ greatly from the KLI 
orbitals. How strong the OEP-SIC localization is will 
depend on the specific system. Generally speaking, we 
expect localization effects to be even more pronounced 
in SIC schemes using orbital dependent potential s 1131 or 
orbital localizing transformations 21, 22, 23, 24], 25|. 

With the ground-state structure of OEP-SIC discussed 
we come to the most important aspect of this manuscript, 



the calculation of the electrical response. As a first test 
we calculated the response of the two dimers Na2 and N 2 , 
which can be seen as representing the "extreme ends" of 
dimer bonding with a soft single and a strong triple bond, 
respectively. The OEP-SIC polarizability (tensor aver- 
age in a§) is obtained as 274 for Na 2 (KLI-SIC performs 
similar) and 10.3 for N2 (no convergence for KLI-SIC). 
The value for the sodium dimer compares favorably with 
the most recent experimental result of 270 [351 ] . the value 
for the nitrogen dimer is too low but not unreasonable 
[3^ . It is a noteworthy observation that OEP-SIC in- 
creases the polarizability (by 12%) for Na 2 (where LDA 
underestimates) while it decreases it (by 18%) for N 2 
(where LDA overestimates), i.e., it works "in the right 
direction" in both systems. OEP-SIC also yields greatly 
improved eigenvalues. For CH4, e.g., OEP-SIC yields 
e HOMO SIC = 14-56 eV, which compares much better with 
the experimental ionization energy of 14.42 eV than the 
LDA value s^oiio =9.52 eV. 

The true and most important test is how OEP-SIC per- 
forms for the response of extended systems where semilo- 
cal functionals fail badly. This is tested by calculating 
the response of the hydrogen chains. The Kohn-Sham 
SIC longitudinal static electric polarizabilities obtained 
from an accurate finite field approach [37I ] are shown in 
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H 8 


H10 


H12 


LSDA 


37.6 


73 


115 


162 


211 


KLI-SIC 


19.4 


60.3 


98.2 


131.7 


193.6 


OEP-EXX 


32.2 


56.6 


84.2 


n/a 


138.1 


OEP-SIC 


30.6 


48.7 


80.1 


98.8 


129.8 


MP4 


29.5 


51.6 


75.9 


n/a 


126.9 



TABLE I: Longitudinal polarizability of hydrogen chains in a§ 
obtained with different exchange-correlation approximations. 
M0ller-Plesset- (MP4) results taken from [38j], exact-exchange 
OEP (OEP-EXX) from 0. KLI polarizabilities were calcu- 
lated from the dipole moment, see discussion in [37| ]. 
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Table [J together with LSDA, exact-exchange OEP (OEP- 
EXX) , and fourth-order M0ller-Plesset perturbation the- 
ory (MP4) results. The MP4 results are close to the exact 
values and serve as the quasi-exact benchmark. The first 
observation is that the KLI-SIC results vary unsystemat- 
ically - the polarizability of H4 is substantially underes- 
timated, whereas the polarizability of all other chains is 
overestimated. Comparison with OEP-EXX and LSDA 
shows that KLI-SIC improves over LSDA, but is less ac- 
curate than exchange-only theory. The picture changes 
when SIC is employed with the true, self-consistent OEP 
instead of with the KLI-approximation: KLI-SIC and 
OEP-SIC polarizabilities are rather different. Compar- 
ing OEP-SIC to the wavefunction based results shows 
that OEP-SIC polarizabilities are within a few percent 
of the MP4 results in all cases and are noticeably closer 
to the MP4 values than the exchange-only OEP results, 
which up to now represented the best density functional 
results for such systems. 

One may wonder why the SIC functional, in which lo- 
calization of the orbitals plays an important role, and 
exact exchange, which is unitarily invariant and thus in- 
dependent of orbital localization, can both lead to a rea- 
sonable description of the chain response. The solution 
lies in the interpretation of the exchange part of the SIC 
functional: The Hartree self-interaction correction corre- 
sponds to the self-exchange part of the EXX functional, 
and it is well known that the diagonal (self-)exchangc 
integrals are the dominant part of exchange, i.e., they 
are noticeably larger than the off-diagonal exchange in- 
tegrals. The larger the diagonal "classical" parts of the 
exchange energy are in comparison with its off-diagonal 
parts, the more accurate becomes the SIC description 
which neglects the off-diagonal parts. Since the diago- 
nal parts are typically maximal for localized orbitals, it 
becomes clear why localized orbitals are crucial in the 
SIC-approach. So from this perspective, SIC takes into 
account the most important part of EXX at the cost of 
needing to employ localized orbitals, but with the huge 
benefit of greatly reducing the number of exchange inte- 
grals that have to be evaluated. In addition, SIC offers 
an improvement over bare EXX that can be attributed 
to the non-EXX parts of the functional. Following [3] 
one can also show that the improved OEP-SIC polar- 
izabilities stem from a field-counteracting-term in the 
response-part of the exchange-correlation-potential [3!| . 
Thus, OEP-SIC is an approach which allows to reliably 
investigate the electrical response of a broad range of 
molecular systems. 

Note added after submission: We learned of studies 
by Pemmaraju, Sanvito, and Burke, and Ruzsinszky, 
Perdew, Csonka, Scuseria, and Vydrov, which also find 
that SIC improves the response. 

S.K. acknowledges financial support by the German- 
Israel Foundation. 



[1] D. J. Tozer, J. Chem. Phys. 119, 12697 (2003). 
[2] N. T. Maitra, J. Chem. Phys. 122, 234104 (2005). 
[3] S. J. A. van Gisbergen et al, Phys. Rev. Lett. 83, 694 
(1999). 

[4] C. Toher et al, Phys. Rev. Lett. 95, 146402 (2005). 
[5] M. van Faassen et al, Phys. Rev. Lett. 88, 186401 (2002). 
[6] N.T. Maitra and M. van Faassen, J. Chem. Phys. 126, 
191106 (2007). 

[7] P. Mori-Sanchez, Q. Wu, and W. Yang, J. Chem. Phys. 

119, 11001 (2003). 
[8] S. Kiimmel, L. Kronik, and J. P. Perdew, Phys. Rev. Lett. 

93, 213002 (2004). 
[9] H. Sekino et al, J. Chem. Phys. 126, 014107 (2007). 
[10] Bulat et al, J. Chem. Phys. 123, 014319 (2005); Cham- 
pagne et al, J. Chem. Phys. 125, 194114 (2006). 
[11] S. Kiimmel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008). 
[12] J.B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A 46, 
5453 (1992). 

[13] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 
(1981). 

[14] P. Cortona, Phys. Rev. A 34, 769 (1986); U. Lundin, and 
O. Eriksson, Int. J. Quantum Chem. 81, 247 (2001); H.- 
J. Unger, Phys. Lett. A 284, 124 (2001); P. Novak et al, 
Phys. Rev. B 67, 140403(R) (2003); O.A. Vydrov et ai, 
J. Chem. Phys. 124, 094108 (2006). 

[15] S. Goedecker and C. J. Umrigar, Phys. Rev. A 55, 1765 
(1997). 

[16] O.A. Vydrov and G. E. Scuseria, J. Chem. Phys. 121, 
8187 (2004). 

[17] J. Chen et al, Phys. Rev. A 54, 3939 (1996). 

[18] CA. Ullrich, P.-G. Reinhard, and E. Suraud, Phys. Rev. 
A 62, 053202 (2000). 

[19] T. Grabo et al, in Strong Coulomb Correlation in Elec- 
tronic Structure: Beyond the Local Density Approxima- 
tion, edited by V. Anisimov (Gordon & Breach, Tokyo, 
2000). 

[20] S.-I Chu, J. Chem. Phys. 123, 062207 (2005). 

[21] A. Svane and O. Gunnarsson, Phys. Rev. Lett. 65, 1148 

(1990); Solid State Comm. 76, 851 (1990). 
[22] W. M. Temmerman et al, Phys. Rev. Lett. 86, 2435 

(2001). 

[23] M. R. Pederson, R. A. Heaton, and C. C. Lin, J. Chem. 
Phys. 80, 1972 (1984); J. Chem. Phys. 82, 2688 (1985); 
M. R. Pederson and C. C. Lin, J. Chem. Phys. 88, 1807 
(1988). 

[24] J. Garza, J. A. Nichols, and D. A. Dixon, J. Chem. Phys. 

112, 7880 (2000). 
[25] S. Patchkovskii, J. Autschbach, and T. Ziegler, J. Chem. 

Phys. 115, 26 (2001). 
[26] J. P. Perdew, Adv. in Quantum Chem. 21, 113 (1990). 
[27] S. Kiimmel and J. P. Perdew, Mol. Phys. 101, 1363 

(2003). 

[28] B. Champagne et al, J. Chem. Phys. 109, 10489 (1998). 
[29] P. Umari et al, Phys. Rev. Lett. 95, 207602 (2005). 
[30] L. Kronik et al, Phys. Stat. Sol. (b) 243, 1063 (2006). 
[31] S. Kiimmel and J. P. Perdew, Phys. Rev. Lett. 90, 043004 

(2003); Phys. Rev. B 68, 035103 (2003). 
[32] M. Griming, O. V. Gritsenko, and E. J. Baerends, J. 

Chem. Phys. 116, 6435 (2002). 
[33] D. Vanderbilt, Phys. Rev. Lett. 79, 3966 (1997); and 

references therein. 



■5 



[34] M. Mundt et al., Phys. Rev. A 75, 050501(R) (2007). 

[35] D. Rayane et al, Eur. Phys. J. D 9, 243 (1999); but 
comparisons have to be done carefully, see S. Kiimmel, 
J. Akola, and M. Manninen, Phys. Rev. Lett. 84, 3827 
(2000); L. Kronik et al, J. Chem. Phys. 115, 4322 (2001). 

[36] The limited performance of SIC for N2 (compared to the 
experimental result of 11.7 given in G. D. Zeiss and W. J. 
Meath, Mol. Phys. 33, 1155 (1977) ) may have its reason 
in SIC taking out too much of the local exchange which 
models non-dynamical correlation. N2 with a triple bond 
may thus represent the "worst case scenario" for SIC. 



For a discussion of SIC and non-dynamical correlation 
see, e.g., V. Polo, E. Kraka, and D. Cremer, Mol. Phys. 
100, 1771 (2002). 
[37] S. Kiimmel and L. Kronik, Comput. Mater. Sci. 35, 321 
(2006). 

[38] B. Champagne et al, Phys. Rev. A 52, 178 (1995); 52, 
1039 (1995). 

[39] T. Korzdorfer, Diploma thesis, University of Bayreuth 
(2006). 



